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NLO QCD corrections to the production of two bottom-antibottom pairs at the LHC 
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We report the results of a computation of the full next-to-leading order QCD corrections to the 
production of two bb pairs at the LHC. This calculation at the parton level provides predictions 
for well separated 6-jets. The results show that the next-to-leading order corrections lead to an 
enhancement of the cross-section for the central scale choice by roughly 50% with respect to the 
leading order result. The theoretical uncertainty estimated by variation of the renormalization and 
factorization scales is strongly reduced by the inclusion of next-to-leading order corrections. 
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INTRODUCTION 



The search for the Higgs boson, and more in general 
the study of the Electroweak Symmetry breaking mech- 
anism, is a major goal of the experiments at the LHC 
collider at CERN. In various extensions of the Standard 
Model the signature of two light Higgs bosons decaying in 
two pairs of 6-quarks, hh — > bbbb, is a viable channel for 
the Higgs Boson discovery. Examples of these models are 
the Minimally Supersymmetric Standard Model (MSSM) 
for large values of tan /3 and moderate ttia [H-Q , hidden 
valley scenarios where the decay of hadrons of an addi- 
tional gauge group can produce additional fr-jets 0, i| 
and two Higgs doublet models. The possibility of mea- 
suring the Higgs self-coupling through H — > hh — > bbbb 
has been investigated in [7[ ■ This and other related stud- 
ies, however show that such a measurement would be 
extremely difficult, primarily due to the large Standard 
Model background. The precise knowledge of the bbbb 
final state within the Standard Model is therefore an im- 
portant factor for the success of these measurements. 

Because of its importance this process has been added 
to the Les Houches wish list of relevant next-to-leading 
order calculations 

In an earlier publication Q , we presented the next-to- 
leading order (NLO) QCD corrections to the production 
of bbbb via quark-antiquark annihilation. In the present 
Letter we complete the existing work including the gluon 
initiated contributions and present the results for the full 
NLO QCD corrections to pp -> bbbb at the LHC. 

We show that the inclusion of the NLO corrections 
reduces the unphysical scale dependence of the leading 
order (LO) prediction greatly, improving the precision 
of this prediction and allowing a better estimation of the 
Standard Model background to possible New Physics sig- 
nals in this channel. 



METHOD 

A complete NLO QCD description requires the calcu- 
lation of the 2 — > 4 subprocesses qq — > bbbb and gg — > bbbb 
at the tree and the one-loop level as well as the 2 — > 5 par- 
ticle processes qq — > bbbbg, gg — > bbbbg and { q ] g — > bbbb ( q } 
at tree level. 

We sum over four massless quark flavours q e 
{u,d,s,c} in the initial state. Neglecting the contribu- 
tion from initial state 6-quarks is justified by the small- 
ness of the b parton distribution function (PDF) with 
respect to the other quark PDFs. Moreover the fact that 
the gluon-gluon channel is the dominant contribution at 
LHC energies further reduces the relative importance of 
the quark channels. We treat the &-quarks as massless, 
which is a very good approximation for LHC kinematics 
also due to the cuts imposed in order for the final state 
6-quarks to be detected and separated in phase space. 
Effects of the heavy top quark arc neglected altogether 
in the final result after having shown that they are nu- 
merically not important. 

The LO and the real radiation matrix elements are gen- 
erated using MadGraph jsj. For the subtraction of the 
infrared singularities we use Catani-Seymour dipoles ficj . 
supplemented with a slicing parameter a as proposed 
in [HI, [Hj], implemented in the MadDipole package 13, 



As described in our earlier work [8j, we compute the 
one loop corrections to scattering matrix elements us- 
ing an approach based on Feynman diagrams. The code 
for the numerical evaluation of the virtual corrections is 
generated using the automated one-loop matrix element 
generator golem- 2.0 [l5l - [l7| which employs QGraf 18 1, 
Form [ljj, the Form library Spinney [2(| and the code 
generator Haggies [21| at intermediate levels of the dia- 
gram and code generation. The reduction and evaluation 
of the loop integrals is performed using the Samurai 22] 
and OneLoop [23( packages respectively. 
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The inte grat ion over phase space is carried out using 
MadEvent 24| and it has been split up in independent 
parts in order to optimize the computational time re- 
quired. 

The first contribution consists of the real emission ma- 
trix element supplemented with the subtraction terms. 
The integration of this contribution over the correspond- 
ing 13-dimensional phase space is one of the main compu- 
tational bottlenecks of such a calculation. This integra- 
tion has been performed using up to 2 • 10 9 phase space 
points. For the qq and gg subprocess the evaluation of 
a single phase space point requires the evaluation of 30 
subtraction terms for each partonic channel and an addi- 
tional 10 subtraction terms are needed for the qg channel. 
This means that a substantial fraction of CPU time is 
spent calculating the dipole contributions. In such a sit- 
uation the use of a value smaller than one for the slicing 
parameter a, as proposed in 111 Il2j. speeds up substan- 
tially the computation by avoiding that each subtraction 
term is evaluated for each phase-space point. Besides the 
reduction of the computational time per point this set- 
ting has a second advantage. If not close to a singularity 
the integrand is given just by the real emission matrix 
element. Close to a singularity, where also subtraction 
terms are calculated, these subtraction terms per defi- 
nition have the same kinematical structure as the real 
emission matrix element which is not necessarily true for 
an arbitrary point in phase space. So for each point the 
integrand is either exactly the real emission matrix el- 
ement or something with the same structure but with 
one singularity subtracted. But as this is an integrand 
where our integration routine is optimized for, choosing 
a value for a smaller than one leads to an improvement 
of the convergence of the integral. In our calculation we 
set a = 0.01. 

The second contribution to the integration combines 
the tree-level contribution and the integrated subtrac- 
tion terms. The virtual matrix element is integrated over 
phase space by reweighting a sample unweighted Born 
level events, as described in [lj|. This leads to a con- 
siderable reduction of the required CPU-time since less 
phase space points have to be evaluated. For the results 
shown below, event samples consisting of 10 — 10 5 un- 
weighted events have been used. The LO event samples 
used for the reweighting have been generated with Mad- 



Event [2J and WHIZARD 25, 26] 



In order to establish the correctness of the results ob- 
tained we have performed a number of non-trivial tests. 
The dipole contributions for single phase space points 
and at the phase space integration level have been com- 
pared with the HELAC code 271 l28| and agreement has 
been established up to double precision accuracy for sin- 
gle phase space points and within integration errors for 
the integrated results. The phase space integration of the 
dipole contributions is validated by checking the indepen- 
dence of the result of the slicing parameter a. Also the 



cancellation of the single and double poles between the 
virtual amplitude and the integrated subtraction terms 
has been verified. Finally, the virtual matrix element 
computation for a single phase space point has been com- 
pared to the result published in (29| . In order to perform 
this comparison the contribution from top quark loops 
has been added, even though it is neglected in the re- 
sults presented in the following section. Our result is in 
agreement with the result of [29( , providing a very strong 
test of our virtual contributions computations. 



RESULTS 

In the following we consider the process pp — > bbbb + 
X at the LHC at a center of mass energy of y/s — 
14TeV. The final state jets are defined by apply- 
ing the fcT-algorithm as explained in [30| with a ra- 
dius in i?-space of 0.8. More precisely, the jet algo- 
rithm requires exactly four &-jets in the final state for 
the event to be accepted. All jets are required to lie 
within a rapidity range of 1 77 ( £> ^ ) | < 2.5 and to have 
a transverse momentum px(bj) > 30GeV. We im- 
pose a separation cut between the jets of AR(bi,bj) — 
yj (0i — (f)j) 2 + {rji — rjj) 2 > 0.8. All results have been 
obtained using the CTEQ6M parton distribution func- 
tions [3lj | with two-loop running of a s both for the LO 
and the NLO cross-section evaluations and a s (Mz) — 
0.118. 

The unphysical renormalization and factorization 
scales are usually chosen to be in the vicinity of the typ- 
ical scale of the process. For processes where heavy par- 
ticles such as top-quarks or lU/Z-bosons are involved the 
masses of these particles provide a natural choice. In 
our case, dealing with massless particles only there is no 
such scale, the only scale involved in this process is the 
PT-cut imposed to define the &-jets. In this respect the 
process considered here is similar to the production of 
four light jets. In [12j it has been shown that the aver- 
age transverse momentum px of a jet is a good choice for 
the production of three jets in hadron-hadron collisions. 
On the basis of the scale choice we made earlier for the 
quark initiated case Q , we define the central scale to be 
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(1) 



which turns out to be of the same order of magnitude as 
the average pr of the jets. 

In Figure [T] we plot the variation of the total cross sec- 
tion for the production of two bottom-antibottom pairs 
at the LHC, with the cuts described previously, when the 
renormalization scale /i r and the factorization fj, p are var- 
ied together, with x defined as the ratio to the the central 
scale, /j, r = fj,F = x ■ fiQ. 
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FIG. 1. Total cross section as a function of the scale n = 
j-i r = fiF — x-fio. Renormalization and factorization scale are 
varied in the same direction. 

If we set the renormalization and factorization scale 
to the value /zo as in Eq. (JTJ we find for the total cross 
section with the cuts described above 

= 140 - 48 ± °- 64 P b - ( 2 ) 
This means that for our preferred choice of scales we 
find that the inclusion of the NLO contribution leads to 
an increase of nearly 50% of the total cross section with 
respect to the LO result of a LO = 94.88 ± 0.14 pb. 
However one observes that, as expected, the dependence 
of the result on the unphysical scales is strongly reduced 
in the NLO result with respect to the leading order one. 




mfb^tg [GeV] 

FIG. 2. Invariant mass distribution of the two 6-jets with 
the highest pr- The black shaded area denotes the tree level 
contribution, the red area denotes the NLO cross-section. The 
error bands for both histograms are determined by a scale 
variation between po/2 and 2/j.o. 

In Figure [5] we plot the invariant mass distribution of 
the two 6-jets with the highest transverse momentum. 



The error bands are obtained by a variation of the scales 
between fio/2 and 2/j, . With respect to the LO result 
one observes a shift of the distribution to lower energies 
due to the inclusion of the radiative corrections. The pr 




FIG. 3. pr distribution of the hardest jet. The error bands 
are defined as in Figure [2] 

distribution of the jet with the highest px is shown in Fig- 
ure [3] Here, the radiative corrections enhance the distri- 
bution at higher momentum. Both Figure [2] and Figure [3] 
show a significant reduction of the error induced by scale 
uncertainties on differential distributions. Moreover, the 
distortion in the shapes of differential distributions when 
going from LO to NLO suggests that the application of 
a global if -factor is not sufficient in order to accurately 
describe the higher-order effects. 

A complete phenomenological study of the production 
of two bottom- antibottom pairs at the LHC, including 
the study of PDF uncertainties and the effects of varying 
the cuts on the final state 6-jets is beyond the scope of 
the present Letter and will be the subject of an upcoming 
publication. 

CONCLUSIONS 

We have calculated the next-to-leading order QCD 
corrections to the production of two bottom-antibottom 
quark pairs at the LHC. This calculation has been imple- 
mented in a highly automated framework for the compu- 
tation of NLO QCD corrections (the golem-2 . frame- 
work) which is based on a Feynman diagrammatic ap- 
proach for the evaluation of virtual corrections imple- 
mented in the Samurai and OneLoop packages, inter- 
faced to the Madgraph/Madevent and MadDipole pro- 
grams for the evaluation of the leading-order and dipole 
subtraction contributions and the phase space integra- 
tion. 

The inclusion of the NLO corrections leads to a sig- 
nificant reduction of the uncertainties due to unphysical 
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scale dependence of the LO result, while enhancing the 
cross section by 50% for our central scale choice. Fur- 
thermore, we have shown that the radiative corrections 
lead to changes in the overall shape of the distributions, 
which cannot be accounted for in a reliable way by a 
simple rescaling of the leading order predictions. 

This reduced theoretical uncertainty improves the 
prospects for the use of the bbbb channel in searches of 
Higgs bosons in various extensions of the the standard 
model like SUSY, two Higgs doublet models or hidden 
valley models. 
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